f 


Reft^taduced 
^  iRe 


ARMED  SERVICES  TECHRICAl  INFORMAnON  AGENCY 
ARLINGTON  HALL  STAHON 
AEINGTON  12,  VIRGINU 


NOTICE:  When  govensaent  or  other  drawings,  speci¬ 
fications  or  other  data  are  used  for  any  purpose 
other  than  In  connection  with  a  definitely  related 
govexnment  procurement  operation,  the  U.  S. 
Goyexnment  thereby  Incurs  no  responsibility,  nor  any 
obligation  idiatsoever;  and  the  fact  that  the  Govern¬ 
ment  nay  have  fonulated,  furnished,  or  In  any  way 
supplied  the  said  drawings,  specifications,  or  other 
data  Is  not  to  be  re^urded  by  Implication  or  other¬ 
wise  as  In  any  manner  licensing  the  holder  or  any 
other  person  of  corporation,  or  conveying  any  rl^ts 
or  pezslsslon  to  manufacture,  use  or  sell  any 
patented  Invention  that  may  In  any  way  be  related 
thereto. 


siRucnm 


A  BASIC  SET  OF  MATHEMATICAL  SUBROUTINES  FOR  LARC 


by 


Naomi  C,  MiUet 


March  1961 


Report  1303 


TABLE  OF  CONTENTS 

Page 

ABSTRACT . 1 

INTRODUCTION . . .  .  1 

SUBROUTINES. .  3 

Square  Root .  3 

Cube  Root . . 8 

Arc  Tangent  I . 14 

Arc  Tangent  D  .  .  .  . . 20 

Natural  Logarithm  . .  24 

Sine-Cosine . 30 

Exponential . 36 

Complex  Addition-Subtraction .  42 

Complex  Multiplication .  44 

Complex  Division  . .  46 

Reciprocal  of  a  Complex  Number . 48 

Absolute  Value  of  a  Complex  Number .  50 

Argument  of  a  Complex  Number.  . .  52 

COMMENTS  ON  THE  USE  OF  THE  ASSEMBLY  AND 

SIMULATOR .  55 

APPENDIX  A.  Optimum  Linear  Approximation  Used 

in  the  Square  Root  Routine.  . . 57 

APPENDIX  B.  Derivation  of  the  Iterative  Procedure 

Used  in  the  Cube  Root  Subroutine .  60 

BIBLIOGRAPHY . 62 

LIST  OF  TABLES 

Table  1  ^  Comparison  of  Computed  Square  Root  Values 

with  True  Values . .  .  5 

Table  2  -  Comparison  of  Computed  Cube  Root  Values 

with  True  Values .  10 

ii 


LIST  OF  TABLES 


Table  3  -  Comparison  of  Computed  Arc  Tangent 

Values  (Livermore)  with  True  Values.  ......  16 

Table  4  Comparison  of  Computed  Arc  Tangent 

Values  (Rem  Rand)  with  True  Values .  21 

Table  5  Comparison  of  Computed  Natural  Logarithm 

Values  with  True  Values . .  .  25 

Table  6  -  Comparison  of  Computed  Sine  Values  with 

True  Values. . 32 

Table  7  -  Comparison  of  Computed  Cosine  Values  with 

True  Values . 32 

Table  8  -  Comparison  of  Computed  Exponential  Values 

with  True  Values.  . . 38 


ACKNOWLEDGEMENT 


The  author  wishes  to  thank  Dr.  John  W.  Wrench,  Jr. 
for  the  computations  of  the  true  values  which  are  in  the 
tables  following  each  routine.  His  comments  were  helpful 
in  successfully  completing  this  report. 


ABSTRACT 

This  report  presents  a  basic  set  of  mathematical 
Subroutines  for  LARC  selected  from  those  programmed 
by  personnel  of  the  University  of  California  Radiation 
Laboratory,  Remington  Rand^and  the  Applied  Mathematics 
Laboratory  of  the  David  Taylor  Model  Basin.  Each  subroutine 
was  analyzed,  then  code-checked  on  UNTVAC  I  usir^  LARC 
assembly  and  simulation  routines  (LISA  and  LIS).  Information 
is  given  for  each  subroutine  in  a  form  for  easy  application^ 
Derivations  of  numerical  methods  used,  tables  of  computed 
values  obtained  from  running  the  subroutinej^and  the  true 
values  are  included. 

INTRODUCTION 

Since  a  great  deal  of  programming  effort  will  be  involved  in 
code -checking  large  scale  problems  for  the  LARC,  the  availability 
of  a  reliable  library  of  mathematical  subroutines  will  be  of 
considerable  help. 

One  way  of  code^checking  the  subroutines  in  advance  has  been 

developed  by  Applied  Mathenuitics  Laboratory  personnel  who  have 

programmed  a  routine,  LIS,  ^  which  will  simulate  LARC  computing 

unit  operations  on  UNIVAC  I.  In  addition,  Applied  Mathematics 

Laboratory  personnel  have  developed  an  assembly  language  and 

2 

programmed  an  assembly  routine,  LISA,  which  will  produce  an 
input  program  for  the  simulator,  LIS,  and  for  the  LARC  Computing 
Unit. 


^References  are  listed  on  page  62. 


A  basic  set  of  subroutines  was  selected  from  those  programmed 
for  LARC  by  the  University  of  California  Radiation  Laboratory, 
Remington  Rand,  and  Applied  Mathematics  Laboratory  personnel. 
They  were  analyzed  and  prepared  for  code^checking  on  UNIVAC  I, 
using  LISA  and  LIS.  The  accuracy  of  each  routine  was  checked 
and  compared  with  existing  tables.  The  result  is  a  set  of  subroutines 
which  are  presented  here  in  a  form  for  easy  reference. 

For  each  routine  the  following  information  is  included:  name, 
purpose,  restrictions,  accuracy,  use,  timing,  programmer, 
method,  cases  used  in  code-checking,  and  coding  in  the  form 
required  as  input  for  LISA.  The  address  of  the  memory  location 
from  which  a  subroutine  is  entered  is  stored  in  register  01.  The 
exit  from  each  subroutine  is  a  return  to  the  address  following  that 
stored  in  register  01. 
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SUBROUTINES 


SQUARE  ROOT 
File  Na^e:  SQR 

Purpose:  To  calculate  the  square  root  of  a  nonnegative,  single- 
precision  floating-point  number  x. 

Accttfacy:  ''HT is  obtained  to  9  significant  figures. 

Usage: 

(a)  Entrance 

1.  Starting  location:  SQR. 

2.  Argument:  ccmtents  of  register  02 » 

(b)  Storage  Requirements 

1.  Five  consecutive  registers:  01  through  05 » 

2.  Instructions:  22  lines  beginning  at  location  SQR# 

3.  Constants:  24  lines  beginning  at  location  CON# 

4.  Error  retwn:  location  SQN. 

(c)  Exit 

Res  alts  stored  in  register  02  . 

(d)  Error  Detected 

If  X  <  0  ccMitrol  is  sent  to  memory  location  SQN. 
Timing:  160yus 

Based  On  Program  By:  L.  Barr,  Livermore 
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Methods  The  determination  of  the  square  root  Of  x  proceeds  as  follows: 

(a)  If  X  =  0,  the  result  is  zero. 

(b)  If  X  <  0,  there  is  an  error  return. 

(c)  If  X  >  0,  then 

(1)  X  =  M  •  10^  where  N  is  an  mteger  and  1  <  M  <  10,  and 
(2r^  is  determined  as  follows. 

The  initial  approximation  to  'yjM  is  obtained  using  a 
set  of  straight  line  segments  to  approximate  y  n^M  <  n+l 

where  n  +  l,  . . . ,  9.  The  equation  for  the  line  Segment  for  the 
nth  interval  is  considered  to  be  written  aS  follows: 
y*  =  a^M+bn 

The  a}^  and  b^  are  obtained  from  stored  tables.  The  determination 
of  and  b^  is  such  that  j  minimized  for  the 

interval  n_<  M  ^n+l.  ♦  This  approximation  is  accurate  to  three 
significant  figures.  Two  Newton  iterations  are  then  used  to 
determine  the  square  root  to  nine  significant  figures.  ^ 

The  square  roots  of  the  numbers  shown  in  Table  1 
have  been  used  in  checking  out  SQR.  In  each  case  nine  significant 
figures  were  obtained. 


♦See  APPENDIX  A. 
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X 


0. 00000000 
1.00000000 
2. 00000000 
3.00000000 
4.00000000 
5. 00000000 
6. 00000000 
7. 00000000 
8. 00000000 
9. 00000000 
26. 5938654 
209. 000000 
0. 00265938684 
0.  000265938654 .Q 
0. 999999999x10 


Computed 


0. 00000000 
1.00000000 
1.41421356 
1. 73205080 
2.00000000 
2. 23606797 
2. 44948974 
2.  64575131 
2. 82842711 
3.00000000 
5. 15692402 
14. 4568322 
0. 0515692402 
0. 0163076256  , 
0. 316227765x10 


True  Value  ofVjT 


0. 00000000 
l.OOOOOOOO 
1.41421356 

1.  73205081 

2.  OOOOOOOO 
2.23606798 
2.44948974 
2.64575131 
2.82842712 
3.00000000 
5. 15692403 

14. 4568323 
0.  0515692403 
0. 0163076256 
0. 316227766x10 


TABLE  1 

Comparison  of  Computed  Square  Root  Values  with  True  Values 
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DUiCB 


LAKC  C.U.  ASSEMBLER 
CODINO  SHEET 


Problem  SQR 


Page:  1 
Date: 

Programmer:  L.  Barr 


LARC  C.U.  ASSEMBLER 
CODING  SHEET 


Problem  SQR 


CUBE  ROOT 


File  Name:  CUR 

Purpose;  To  calculate  the  cube  root  of  a  single^precision, 
floating-point  humbe^',  X. 

Accuracy:  ^  X  is  obtained  to  9  significant  figures. 


(a)  Entrance 

1.  Starting  location:  CUR. 

2.  Argument:  contents  of  register  04. 

(b)  Storage  requirements 

1.  Four  consecutive  registers:  01  through  04, 

2.  instructions:  48  lines  beginning  at  location  CUR  > 

3.  Constants:  10  Imes  beginning  at  location  CQN. 

4.  Working  storage:  5  lines  beginning  at  location  COM 

(c)  Exit 

1.  Results  stored  in  03. 

2.  Argument  stored  in  02. 

(d)  Errors  detected 

None . 

Timing:  SOOyUs 

Based  on  Program  By:  Robert  Shafer,  Livermore 
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Mettipd:  The  determination  of  the  cube  root  of  X  proceeds  as  follows: 

(a)  If  X  =  0^  then  X^/^  is  zero. 

(b)  H  X  0,  then  X  =  M  •  10^  where  N  is  an  integer  and  .  1  <  M  <  1. 0 
Therefore, 


which  is  computed  as  follows. 

N 

(1)  The  division  of  ^  is  performed.  Then  we  have 
3 

where  I  is  an  integer  and  J  =  0,  1/3  or  2/3,  and 


loN/3  ^  iqI  .  iqJ 


^0 


,  and  10^/^  ar  e 


is  computed.  The  numbers  10  , 
stored  as  constants  in  the  Subroutine. 

(2)  To  compute  three  iterations  are  performed; 


using  the  formula 

™(n+l)  "  ( 


3 

2M  +  mn 
M  +  2mn^ 


The  initial  approximation  mQ  =  .  7  is  used  to  obtain  m^.  The 

1/3 

third  iteration  yields  an  approximation  to  M  ^  which  is 


accurate  to  the  prescribed  number  of  significant  figures.  A 
derivation  of  this  iteration  procedure  is  given  in  Appendix  B. 
The  results  of  code -checking  may  be  found  in  Table  2. 
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X 

o 

Computed 

:  True  Value  of 

0. oooooooooooo 

;  0. 00000000 

1  0.00000000 

1.000000000000 

1.00000000 

1.00000000 

i. 500000000000 

1  1.14471424 

1.14471424 

2. oooooooooooo 

1.25992104 

1.25992105 

2. 500000000000 

1.35720880 

1.35720880 

3.000000000000  i 

1  1.44224957 

1.44224957 

3.500000000000 

;  1.51829448 

1.51829448 

4.  oooooooooooo 

1.58740105 

1.58740105 

5.000000000000 

1.70997594 

1.70997595 

6.000000000000 

1.81712059 

1.81712059 

7. oooooooooooo 

!  1. 19293118 

1.91293118 

-0.000367905300 

1-0. 0716547973 

-0. 0716548099 

791. OOOOOOOOOOOO .Q 
0. 999999999x10^^ 

9.24821859 

0. 215442679x10^ 

.  9.24823438 
0.215443469x10^^ 

TABLE  2 

Comparison  of  Computed  Cube  Root  Values  with  True  Values 

< 
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Problem 


LAAC  aU.  ASSEMBLER 
CODING  SHEET 


Page:  2 
Date: 
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DTMB 

LARC  C.U.  ASSEMBLER 

PiBge:  3 

COmNG  SHEET 

Date: 

Programmer: 

Problem 

AnC  TANGENT  I 

File  Name:  AkT 

Purpose:  To  calculate  the  arc  t^gent  of  a  single^precision, 
floating-point  number,  X. 
nestrictions: 

10'^®  <  xj<  10^^ 

Accuracy:  Arc  Tan  X  is  obtained  to  9  significant  figures. 

Usage: 

(a)  Entrance 

1.  Starting  Location:  ART. 

2.  Argument:  contents  of  register  04. 

(b)  Storage  requirements 

1.  Four  consecutive  registers:  01  through  04. 

2.  Instructions:  44  lines  beginning  at  location  ART. 

3.  Constants:  20  lines  beginning  at  location  CON. 

4.  Working  storage:  4  lines  beginning  at  location  COM. 

(c)  Exit 

Results  stored  in  register  02. 

(d)  Errors  detected 
None. 

Timing:  288^ s 

Based  on  Program  By:  Robert  Shafer,  Livermore 
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Method: 


Case  I:  X  >  1 
In  this  case 

tan"^  X  =  tan"^  (1)  +  tan"^^  +  tan"^  1/ 
is  used  where  ^  some  integer  N,  0  <  H  <  10  andt^<  .  1. 

The  details  follow. 


First 

taii"^  X  =  tan"^  (1)  +  tan"^  (  ^^-7-) 
is  used  and  then 
.-1  . 


tan' 


Where 


P  *  )  =  +tan"^ 

X+1 


/  X-1  V 

( - )  +  .  05 

X+1 


the  largest  integer  <  ^  • ) 


Finally  letting  - 
Case  II: 


i 


A 


o(2+^/3>i 


designates 


=  2y,  then  y  <  .1. 


X  =  1 

tan"^  (1)  is  available  directly. 

Case  ni:  0  <  X  <  1 
In  this  case 

tan'^  X  =  tan'^  c?»,  +  tan”^  V 

N 

is  used  where  ai  =  — —  for  some  integer  N, 

10 

0  <  10  andi/<  .  1. 
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Case  IV:  X  <  0. 

In  this  case 

tan"^  X  =  “tan  ^  y,  where  y  =  -X,  y  >  0  is 
used.  Then  Case  I,  H^or  HI,  whiGhever  is 
appropriate,  is  used  to  deterniine  tan”^  y. 

The  arguments  used  in  eode-Ghecking  are  listed  in  Table  3 . 


X 

Arc  Tan  X 

True  Value  of 
Arc  Tan  X 

2.50000000 

1. 19028994 

1. 19028995 

-3.50000000 

-1.29249666 

-1.29249667 

4. 50000000 

1.35212738 

1.35212738 

5.50000000 

1.39094282 

1.39094283 

6.00000000 

1.40564764 

1.40564765 

7. 00000000 

:  1.42889927 

1.42889927 

-0.11000000 

“0. 1095595267 

-0. 1095595268 

0.577x1014 

1.57079632 

1.57079633 

TABLE  3 


Comparison  of  Computed  Arc  Tangent  Values  (Livermore) 

with  True  Values 
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Date; 
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File  Name:  ARC 

Purpose:  To  calculate  the  arc  tangefit  of  a  sihgie^precision,  floatings 
point  nuMber, 

Accura.cy:  Arc  tan  X  is  obtained  to  9  significant  figures. 


(a)  Entrance 

1.  Starting  location:  ARC. 

2.  Argument:  contents  of  register  06. 

(b)  Storage  requirements 

1.  Six  consecutive  registers:  01  through  06. 

2.  Instructions:  37  lines  beginning  at  location  ARC. 

3.  Constants:  10  lines  beginning  at  location  CON. 

(c)  Exit 

ilesults  stored  in  06. 

(d)  Errors  detected 

None. 

Timing:  144yas 

Based  on  Program  By:  A.  B.  Tonik,  Remington  Rand 
Method: 


Case  I:  0  <  X  <  1 


hi  this  case 


arc  tan  X  =  g  +  arc  tan  Z 


where  Z  = 


X  -(>f2  -  1) 

1  +  (V^-  1)  X 
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case  11:  |X|i>l 
In  this  case 

3  It 

arc  tan  X  =  -  ^  -  +  arc  tan  Z, 

8 

where  Z  = 

X  (  fT-  1)  -  1 

The  arguments  used  in  code-checking  are  listed  in  Table  4. 


X 

Arc  Tan  X 

True  Value  of 
Arc  Tan  X 

2.  5 
-3.5 

4.5 

5.5 

6.0 

7.0 

0.  577xl0j^ 
-0.577x10^^ 

1.19028994 

-1.29249666 

1.35212738 

1.39094282 

1  1.40564129 
1.42888707 
-0.0996686525 

1. 50920769 
-1.50920769 

1. 19028995 
^  -1.29249667 

1.  35212738 

1.  39094283 

1. 40564765 

1.  42889927 
-0.0996686525 
1.57079633 
-1.57079633 

TABLE  4 

Comparison  of  Computed  Are  Tangent  Values  (Rem.  Rand) 

with  True  Values 
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NATUKAL  LOGAHITHM 


File  Name:  LOG 

Purpose:  To  compute  the  natural  logarithm  of  a  floatihg-^point, 
single-precision  number ,  X. 

Restrictions:  X  >  0 

Accuracy:  In  Xis  obtained  to  9  significant  figures. 

Usage: 

(a)  Entrance 

1.  Starting  location:  LOG. 

2.  Argument:  contents  of  r  egister  02 . 

(b)  Storage  Requirements 

1.  Five  consecutive  registers:  01  through  05. 

2.  histructions:  25  lines  beginning  at  location  LOG. 

3.  Constants:  76  lines  beginning  at  location  LGC. 

(c)  Exit 

Results  stored  in  register  02. 

(d)  Errors  Detected 

If  X  <  0,  control  is  sent  to  memory  location  NLA. 
Timing:  148  yUs 

Based  On  Program  By:  Leota  Barr,  Livermore 
Method:  The  determination  of  In  x  proceeds  as  follows: 

(a)  If  X  <_0,  there  is  an  error  return 

(b)  If  X  >  0, 
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then 


In  X  =  I  In  10  +  In  F 


where  I  is  an  integer  and  ,  1  <  F  <  1. 0 


Let 


where  d  is  a  set  of  stored  constants  whose  use  is  determined 


by  the  first  digit  of  F  such  that  .  5  <  Fd  <  1  #  0 


Let 

F' 


G 

H  +  1 


where  G  is  the  two  most  significant  digits  of  F*  0. 00999999999. 
Then  H  =  <  .  02  and 

p, 

In  (H+l)  4  Cj  H  +  Cg  4  Cg 

The  C|'s  are  derived  using  Chebyshev  polynomials. 

The  arguments  used  in  code^checking  this  subroutine  are  shown  in 


X 

LN  X 

True  Value  of 
LN  X 

0. 001 

-6.90775526 

-6.90775528 

1.5 

0.405465108 

0.405465108 

2. 0 

0. 693147189 

0.693147181 

3.0 

1. 09861230 

1.09861229 

4.0 

1. 38629437 

1.  38629436 

5.0 

1. 60943793 

1.60943791 

6.0 

1.79175947 

1. 79175947 

7.0 

1.94591015 

1.94591015 

8.0 

2. 07944156 

2.07944154 

9.0 

2. 19722459 

2. 19722458 

112.0 

4.71849888 

4.71849887 

40000. 0 

10. 5966348 

10. 5966347 

TABLE  5 


Comparison  of  Computed  Natural  Logarithm  Values  with  True  Values 
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SINE -COSINE 


File  Name:  SINCOS 


Purpose:  To  compute  the  sine  or  cosine  of  a  single-precision,  floating¬ 
point  number,  X,  in  radimis. 

Restrictions:  |x|<10^® 

Accuracy:  Sine  or  cosine  is  obtained  to  9  significant  figures. 


Usage: 

(a)  Entrance 

1.  Starting  location:  SIN  (or  COS). 

2.  Argument:  contents  of  register  02, 

(b)  Storage  requirements 

1.  Three  consecutive  registers:  01  through  03. 

2.  Instructions:  69  lines  beginning  at  location  COS  or  SIN- 

3.  Constants:  22  lines  beginning  at  location  CON. 

4.  Working  storage:  1  line  at  location  COM. 

(c)  EaU 

Results  of  SIN  or  COS  stored  in  03. 

(d)  Errors  detected 

None. 

Timing:  348 yus  average 

Based  On  Program  By:  K.  B.  Williams,  Livermore 


30 


Method: 

Case  A  :  Cosine  Routine 

1  quadrant  storage  location,  QSL 
Case  B:  Sine  Routine 

0  quadrant  storage  location,  QSL 

JT 

Case  I:  0  <  x  X 

In  Cases  I  and  IV  if  x  <  0. 001,  then  x  •  sin  x. 

If  X  >  0. 001  then  sine  or  Cosine  series  is  used,  depending  on 
the  contents  of  location  QSL.  The  coefficients  in  the  series 
are  obtained  from  a. Los  Alamos  SGientific  Laboratory 

4 

report. 

Case  II ;  a  >  0;  x  >  -|- 
Then 

X  •  -f:_  s  Y,  where  Y  consists  of  an  integral  part,  N, 
n 

and  a  fractional  part. 

Then 

N  ^  (QSL)  ^QSL, 

(Y  -  N) 

it 

and  Case  I  is  used. 

„  n 

Case  III;  x  <  0;  ^ 

2  +  QSL  -^QSL 
then  Case  I  is  used. 

Ca.f  IV;  X  <  0;  I  x<>  y 

Cast'  HI  and  Case  11  are  used. 

% 
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Arguments  used  for  checking  out  the  routine  are  shown  in 


Table  6  and  Table  7. 


X 

Sine  X 

True  values 
of  Sine  x 

0.  7 

0. 644217693 

,  0.644217687 

-0.51196 

-0. 489886890 

!  -0. 489886887 

1.7482 

0.984305147 

>  0.984305196 

0.0009 

0. 000900000 

0. 000900000 

i  4. 18879 

-0. 866025309 

-0. 866025301 

5.75959  ' 

-0.  499997010 

-0. 499996996 

-1.05883 

-0.871782915 

-0.871782905  1 

-2. 01586 

-0.902583260 

-0.902583255 

-4. 71239 

1.000000000 

1.000000000 

-5.58505 

0. 642790380 

0. 642790372 

Comparison  of  Computed  Sine  Values  with  True  Values 


X 

Cosine  x 

True  Values 
of  Cosine  x 

0.7 

-0.51196 

1.7482 

0.  0009 

4. 18879 
5.75959 
-1.05883 
-2.01586 
-4.71239 
-5.58505 

0. 764842190 

0. 871786010 
-0. 176474592 

0. 999999600 
-0. 500000190 
0.866027140 
0.489892410 
-0.430515351 

0. 101685500x10  - 
0. 766042130 

0. 764842187 
0.871786004 
-0. 176474593 
0.999999595 
-0.500000177 

0. 866027138 

0. 489892403 
-0.430515352 

0. 101961531xl0‘^ 

0.  766042125 

TABLE  7 

Comparison  of  Computed  Cosine  Values  with  True  Values 
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LARC  C.U.  SIMULATOR 


CODDfG  SHEET 


Problem 


EXPONENTIAL 


File  Name;  EXP 

Purpose;  TO  calculate  the  e;..  onential  e^  of  a  floating-point^  single- 
precision  number,  x. 

Restrictions:  -111  <  x  <  112 

Accuracy:  e^  is  obtained  to  d  significant  figures  . 

Usage: 

(a)  Entrance 

1.  Starting  location:  EXP. 

2.  Argument:  ccmtents  of  register  02. 

(b)  Storage  Requirements 

1.  Eight  consecutive  registers:  01  through  08. 

2.  Instructions:  28  lines  beginning  at  location  EXP. 

3.  Constants:  38  lines  beginning  at  location  CON. 

(c)  Exit 

Results  stored  in  02, 

(d)  Errors  Detected 

1.  If  X  <  -111.0,  then  an  absolute  zero  will  be  stored  in  02. 

2.  If  X  >  112.0,  then  a  transfer  to  memory  location  ETL 
will  be  made. 

Timing:  208yUs 

Based  On  Prograni  By:  Robert  Shafer,  Livermore 
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Method;  The  determination  of  the  exponential  function  proceeds  as 
follows: 

(a)  If  X  ^-111. 0,  then  e*  =  0. 

(ib)  If  X  ^  112. 0,  there  is  an  error  return. 

(g)  If  -111.  0  <  X  <  112. 6,  the  following  two  cases  must 
be  Gonsidered: 

(1)  If  X  <  0 
e*  -  e*  .  e^ 

where  I  is  an  integer  and 
.  1  <  F  <  LO 


Now 

X  =  (I  In  10  +  F  In  10) 
Let 


1 

In  10 


r.  F’ 

^  ‘irni 


F'  =  FlnlO 
e  F'  = 

where  G  =  first  two  digits  of  F' . 


Let  H  =  F'  -  G 


Then 

e®  ~  Cq  +  ?!  H  +  C2  +  C3  +  C4 

5 

where  the  cj's  are  derived  using  Chebyshev  polynomials. 
(2)  If  X  >  0 
then 

gX  _  gI+1  e  F-1 
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Again  let  F  =  ^  ,  andi  the  procedure  is  followed  as 

In  10 


in  ( 1)  above . 

the  arguments  used  in  code -checking  are  shown  in  table  8. 


X 

true  Value  of  e* 

1.0 

2.71828181 

2.71828183 

7.  38905612 

7.38905610 

20. 0855369 

20. 0855369 

4. 0 

54.5981495 

54. 5981500 

5.0 

148.413160 

148.413159 

6.0  1 

403. 428795 

403.428793 

7.0 

1096.63316 

1096.63316 

8.0 

2980. 95795 

2980. 95799 

9.0 

8103.08379 

8103.08393  .. 

111.0 

0. 160948712x10^® 

0. 160948707x10^® 

-110.0 

0. 168891188x10"^^ 

-  ^ 

0. 16889 1188x10“^^ 

tABLE  8 

Comparison  of  Computed  Exponential  Values  with  true  Values 
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COMPLEX  ADDITION-SUBTkACTION 

File  Name:  ADC  SBC 

Purpose;  To  add  or  subtract  two  floatingHpoint,  single^recisionv 
complex  numbers  Zj  =  Xj  +  1X2  and  =  Vj  +  iy2; 
Accuracy:  The  result  of  the  ad<^tion  or  subtraction  is  obtained  to 
9  significant  figiures. 


(a)  Entrance 

1.  Starting  location:  ADC  (or  SBC). 

2.  Each  complex  number  is  represented  in  two  adjacent 
registers,  02  and  03,  04  and  05  containing  the  real 
part  and  the  imaginary  part,  respectively. 

(b)  Storage  Requirements 

1.  Five  consecutive  registers:  01  through  05. 

2.  Instructions:  6  lines  beginnhig  at  location  ADC  (or  SBC). 

(c)  Exit 

1.  Real  psu^t  stored  in  register  02. 

2.  Imaginary  part  stored  in  register  03. 

(d)  Errors  Detected 

None. 

Timing:  16/(s 

Programmed  by:  Naomi  Millet,  David  Taylor  Model  Basin 


Method: 


DTMB 


LAAC  G.U.  ASSEMBLER 
CODING  SHEET 

Problem  ADC 


Page:  1 
Date: 

Programmer:  N.  Millet 
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COMPLEX  MULTIPLICATION 


File  Name:  CPM 

Purpose;  To  multiply  two  floating-point,  single-preeision  ,  complex 
numbers,  Zj  =  +  iXg  and  =  Yj  +  iy^- 

Accuracy  ;  The  result  of  the  multiplication  is  obtained  to 
9  significant  figures. 

Usage; 

(a)  Entrance 

1.  Starting  Location;  CPM. 

2.  Each  complex  number  is  represented  in  two  adjacent 
registers  02  and  03  containing  the  real  part  and  the 
imaginary  part,  respectively. 

(b)  Storage  Requirements 

1.  Five  consecutive  registers;  01  through  05. 

2.  Instructions;  8  lines  beginning  at  location  CPM. 

(c)  Exit 

1.  Real  part  stored  in  register  02. 

2.  Imaginary  part  stored  in  register  03. 

(d)  Errors  Detected 

None. 

Timing;  56^4<s 

Programmed  By;  Naomi  Millet,  David  Taylor  Model  Basin 
Method; 

(Xj  +  iX^)  (  Yj  +  iy^j)  =  (Xj  Yi  -  X2  y^)  +  i(X2  Yi  +  Xj  y2) 
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LARC  C.U.  ASSEMBLER 
CODING  SHEET 


Page:  1 
Date: 

Programmer:  N.  Millet 


Problem  CPM 
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COMPLEX  DIVISION 
File  Name:  GPD 

purpose:  To  ^vide  two  floatingpoint^  singlepf ecision,  complex 
numbers  =  X|  +  iX^  and  =  Vji  +  iy^- 
Accuracy;  The  result  of  the  (^vision  is  obtained  to  9  significant  figiu'es. 
Usage: 

(a)  Entr  ance 

1.  Starting  location:  CPD. 

2.  Each  complex  number  is  represented  in  two  adjacent 
registers,  02  and  03,  04  and  05,  containing  the  real 
part  and  the  imaginary  part,  respectively. 

(b)  Storage  Requirements 

1.  Five  consecutive  registers:  01  through  05. 

2.  Instructions:  16  lines  beginning  at  location  CPM. 

3.  Working  storage:  4  lines  beginning  at  location  COM. 

(c)  Exit 

1,  Real  part  stored  in  02. 

2.  Imaginary  part  stored  in  03. 

(d)  Errors  Detected 

None. 

Timing:  148^s 

Programmed  By:  Naomi  Millet,  David  Taylor  Model  Basin 
Method: 

(X^  -K  iXg)  Xiyi  +X2  y2 
(yi  +  iy2)  Vl^  +  yz^ 
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X2  yi  '  Xi  y2 

y +  y2^ 


DTMB 


LARC  C.U.  ASSEMBLER  Page:  1 

CODIKG  SHEET  Date: 

Programmef: 

Problem  CPD 
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nECIPKCX:AL  OF  A  COMPI^X  NUMBBK 


File  Name;  ftCG 

Piurpose:  Td  calculate  the  reciprocal  of  a  lloating-point.  sihgle- 

precision,  Gorhplex  number,  Z  =  X  +  iX^. 

1  ^ 

Accuracy:  The  reciprocal  is  obtained  to  9  significant  figures. 
Us^e: 

(a)  Entrance 

1.  Starting  location:  RGC. 

2.  The  complex  number  is  represented  in  two  adjacent 
registers,  02  and  03  containing  the  real  part  and  the 
imaginary  part,  respectively. 

(b)  Storage  Requirements: 

1.  Four  consecutive  registers:  Ol- through  04. 

2.  Instructions:  10  lines  beginning  at  location  RGG  . 

3.  Working  storage:  2  lines  beginning  at  location  GOM. 

(c)  Exit 

1.  Real  part  stored  in  02. 

2.  Imaginary  part  stored  in  03, 

(d)  Errors  Detected 

None. 


Timing:  104^s 

Programmed  by:  Naomi  Millet,  David  Taylor  Model  Basin 


Method: 

1 

(Xi+  iX2  ) 


+  i 


/  -  X2 

+X2 


2 


) 
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LARG  C.U. 
CODING 


Problem  RGC 


Date: 

programmer:  N.  Millet 
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ABSOLUTE  VALUE  OF  A  COMPLEX  NUMBER 


File  Name;  MGG 

Purpose:  To  calculate  the  modulus  of  a  floating-point,  single - 
precision,  complex  number  Z  =  +  iX^. 

Aecuracy:  |  X  ^  +  iXg  |  is  obtained  to  9  significant  figures. 

Usage: 

(a)  Entrance 

1.  Starting  location:  MGG^ 

2.  The  complex  number  is  represented  in  two  registers, 
02  and  03,  containing  the  real  part  and  the  imaginary 
part,  respectively. 

(b)  Storage  Requirements 

1.  Five  consecutive  registers:  01  through  05. 

2.  Instructions:  28  lines  beginning  at  location  MGG. 

3.  Constants:  24  lines  beginning  at  location  SQG. 

4.  Working  storage:  1  line  at  STO. 

(c)  Exit 

Results  stored  in  register  02. 

Timing:  192  /Us 

Programmed  By:  Naomi  Millet,  David  Taylor  Model  Basin 
Method: 

I  Xi  .  iXj  1= >1x^7^ 
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DTMB 


LARC  C.U.  ASSEMBLER 
CODING  SHEET 


Problem  MGC 


Page:  1 
Date: 

Programmer:  N.  Millet 
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ARGUMENT  OF  A  COMPLEX  NUMBEm 


File  N^e:  PHA 

Purpose;  To  calculate  the  phase  angle  of  a  floating-point,  single- 

precision,  complex  number,  Z  =  X^  +  1X2. 

Restrictions:  10  t  <  10^^ 

Xi  I 

Accuracy:  The  phase  angle  in  radians  is  obtained  to  9  significant 
figures. 

Usi^e: 

(a)  Entr  ance 

1.  Starting  location:  PHA. 

2.  Each  complex  number  is  represented  in  two  adjacent 
registers,  Xi  in  02  and  X2  in  03  containing  the  real 
part  and  the  imaginary  part,  respectively. 

(b)  Storage  Requirements 

1.  Five  consecutive  registers:  01  through  05. 

2.  Instructions:  47  lines  beginning  at  location  PHA. 

3.  Constants:  20  lines  beginning  at  location  ARC. 

4.  Working  storage:  1  line  at  STO. 

(c)  Exit 

Result  stored  in  register  02. 

(d)  Errors  Detected 

None. 

Timing:  336^  s 
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Programmed  By:  Naomi  Millet,  David  Taylor  Model  Basin 
Method: 

0  =  arc  tan 

Xi 

Use  Livermore  arc  tangent  routine. 
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CODMG  SHEET  Date: 

Progranlmer:N 

Problem  PHA 


1 

I 

1 

II 

1 

M 

I 

H 

II 

1 

I  /  / 

PHA 

i 

1 

7  /  / 

S_- 

T  B 

F 

X 

1 

X 

X 

9 

S 

A 

S 

A  APR  A  A 
9,9  +|0  0  3 

t!0  _  !  _ 

Rl'T  • 

1  1  •  ■ 

TO 

lA 

|A  A 

A 

A 

, 

•  •  • 

' '  •  •  • 

1 

0 

• 

• 

0  0  +  0  0  1 

,  ■  ■  ■ 

i 

1 

1 

1  * 

0  1 

.  •  • 

•  • 

; _ ii_J 

• 

• 

-  1 

' 

• 

_ 1 

• 

1 

1 

1 

1 

. 

1 

1 

' 

1 

1 

1 

_ J _ 

■ 

' 

1 

1 

1 

1 

1 

H 

Millet 


54 


COMMENTS  ON  THE  USE  OF  THE  ASSEMBLY  AND  SIMULATOR 

For  the  types  of  routines  included  in  this  report,  the  LARC  Input 
Simulator  Assembly  (LISA)  assembled  input  for  the  LARC  Instruction 
Simulator  (LIS)  and  providedi  an  analyzer  which  listed  programming 
errors  in  approximately  fifteen  minutes  of  UNIVAC  I  running  time. 
Multiply^defined  symbols  were  not  indicated!,  and  any  new  definitions 
of  a  symbol  previously  defined  were  ignored.  Memory  space  was 
reserved  for  all  registers  from  register  zero  up  to  and  including  the 
highest-numbered  register  mentioned  in  the  specific  program.  Register 
zero  always  contains  zeros.  Other  registers  cannot  be  assumed  to 
contain  zeros  initially. 

LIS  consists  of  four  parts.  Part  1  selects  the  Order-Subroutines 
which  are  required  to  simulate  the  given  routine.  These  order-subroutines 
are  in  x-1  assembly  language.  Part  2  conveHs  these  chosen  order - 
subroutines  into  UNIVAC  machine  language.  Part  3  performs  the  actual 
simulation  of  the  LARC  instructions.  Part  4  edits  the  results  and 
prepares  the  output.  The  running  of  Parts  3  and  4  of  LIS  took  approximately 
five  minutes  for  each  of  the  routines  included  in  this  report.  The  accuracy 
of  the  computed  answers,  permitting  up  to  nine  significant  figures  in  single¬ 
precision  floating-point  and  eleven  significant  figures  in  single-precision 
fixed-point,  may  be  judged  by  the  table  following  each  routine. 
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LIS  provides  the  programmer  with  a  tape  containing  the  input 
to  Part  3  of  LIS.  Therefore,  most  corrections  necessary  during 
code -checking  may  be  made  on  this  tape,  eliminating  the  need  to 
assemble  and  perform  Parts  1  and  2  of  LIS  more  than  once.  However, 
care  should  be  taken  in  correcting  and  adding  to  this  tape  since  it 
will  eontain  only  the  order -subroutines  required  for  simulating 
those  instructions  which  are  in  the  original  input  to  Part  1.  A 
chart  is  contained  in  the  report  on  LIS  listing  the  set  of  instructions 
which  each  order  ^subroutine  can  simulate.  This  will  allow  the 
programmer  to  determine  which  order -subroutines  have  been 
included  in  the  tape  containing  the  input  to  Part  3. 

At  the  writing  of  this  report,  the  OCX  instruction  could 
not  yet  be  handled  by  LIS,  but  it  is  expected  that  this  instruction 
will  be  added  in  the  near  future. 
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APPENDIX  A 


Optimurn  Linear  Approximation  Used  in  the  Square  i.^oot  iloutine 


The  probiem  is  to  approximate  the  curve  y  =  £(x)  by  a 

straight  line  y*  s  ax  +  b  in  an  interval  x  <  x  ^Xj,  x  >  Q, 

over  whlGh  f  ’  (x)  >  0  ^Jid  £”  (x)  <  0. 

(I)  To  determine  the  line  y**  =  a*x  +  b*  throiigh  the  points 

(xq,  ftxg)),  (xj,  f(Xj)),  we  have 

=  X  ^  XQ  f(xi)  -  XI  f(xQi) 

f (Xg)  +  KXj)  *0  '  XJ 


where 


b*  =  ^  ^(^l)  -  H  ^(xq) 

Xq  -  XI 

(2)  To  determine  the  maximum  value  of  (y  -  y**)  in  the 

interval  Xq  <  x  ^Xj^,  x  ^0,let 

d(y-y**)  =  0 

dx 

Since 

y  -  y**  =  f(x)  -  a*x  -  b* 

then 

g  fciL"-)  =  V  (x)  -  a‘. 
dx 

The  maximum  value  is  attained  at  the  point  x*  for  which 
V  (x*)  =  a* 
or 


X*  ^  (a*) 
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max  (y  -  y**)  =  f(x*)  -  (a*  x  ♦  +  b*)  =  d* 

(3)  By  choosing  the  straight  line  parallel  to  y  ** 
y*  =  ax  +  b 
for  which 

max  (y-y*)  =  ^  d* 

where 

a  =  a* 

1 

b  =  b*  +  -  d* 

b*  a*  X*  f  (x*) 

=  -  (  b*  -  a*  X*  +  f(x*)* 

Under  these  circumstances 

(Y-y»)  =  d^  at  X  =  x.,  Xj;  (y-y»)  =  -  d*  at  x  =  x*, 
and  we  have  the  optimum  straight  line  to  approximate  f(x)  in 
the  specified  interval. 

Now 

y  = 

1 

a  -  - - - 

b  ^  7  (  ^  ^1  '  "^0  -  a*  X*  +  f(x*)) 

^0  * 

1  /  ^0  ~  .  X  + 

2  '  Xq  -  Xj  4a*  2a* 
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OPTIMUM  LINEAR 


n 

a 

b 

max 

1 

.4142 

.  5947 

.  009 

2 

.3178 

.  7826 

.004 

S 

.  2679 

.  8306 

.003 

4 

.  2361 

1.0574 

.002 

5 

.2134 

1. 1702 

.001 

6 

.  1963 

1.2729 

.001 

7 

.  1827 

1.3678 

.001 

8 

.  1716 

1.4565 

.001 

9 

.  1623 

1.5400 

.001 

4  (Xq  -  xj) 


UGRL 


n  a  b  max  y^y* 


1 

.  4131 

.  596 

.  009 

2 

.3175 

.784 

.004 

3 

.  2678 

.931 

.  003 

4 

.2358 

1.059 

.002 

5 

.  2133 

1. 170 

.002 

6 

.  1962 

1.273 

.001 

7 

.  1827 

1.368 

.001 

8 

.  1716 

1.457 

.001 

9 

.  1623 

1.540 

.001 
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To  eliminate  x  from  the  right-hand  side  of  this  equation,  we  use  the 
first-order  approximation  obtained  by  terminating  the  Taylor  expansion 
of  f(x)  after  the  second  term,  namely 


Substituting  this  expression  in  Equation  (6)  we  have 


(Xq) 


or 

(9)  X  =Xo  -  . 

2[l'(iio)J  ^  -  t  (V  f"  (xo) 

Note  that  if  f"(xQ)  =  0,  Equation  (9)  reduces  to  the  Newton -Jlaphson 


formula 


X  =  X 


0 


- ' 


f  (Xj,) 


Since  from  Equation  (1)  we  have  f  '  (x)  =  3x*  and  f”  (x)  =  6x, 


we  can  substitute  these  expressions  in  Equation  (9),  obtaining 
X  =  xo  - 


2  (SxJ  )  -  (x®  -  N) 


which  reduces  to 


2N  +  xt 

X  =  ■  V  .  X 

2  xg  +  N  C 


This  is  the  proeedure  used  in  the  Cube  Root  iaoutine,  with 
Xq  =.7  in  the  initial  approxiraation. 
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